#!/usr/bin/python
import pylab
import numpy
import scipy
import os
import cmath
import math


#digit revert order
def dro(index, N):
    width = math.log(N, 2)
    mid = (width-1)/2.0
    n = 0
    for j in range(int(width)):
        bit = (index & (1 << j)) 
        if j<width/2:
            bit = bit << int(((mid-j)*2))
        else:
            bit = bit >> int(((j-mid)*2))
        n = n | bit   
        #print '\t', n
    return n


def bf1(input):
    length = len(input)
    x=[0]*length
    #print 'bf1'
    for i in range(length/2):
        x[i]          = input[i] + input[i+length/2]
        #print '%d=%d+%d'%(i,i,i+length/2)
        x[i+length/2] = input[i] - input[i+length/2]
        #print '%d=%d-%d'%(i+length/2,i,i+length/2)
    return x

def tf1(input):
    length = len(input)
    x=input[:]
    for i in range(length*3/4, length):
        x[i] = ((-1)*1j)*input[i]
        #print '%d multiplied -j'%(i)
    return x

def bf2(input):
    length = len(input)
    x=[0]*length
    #print 'bf2'
    for i in range(length/4):
        x[i]                = input[i]           + input[i+length/4]
        #print '%d=%d+%d'%(i,i,i+length/4)
        x[i+length/4]       = input[i]           - input[i+length/4] 
        #print '%d=%d-%d'%(i+length/4,i,i+length/4)
        x[i+length/2]       = input[i+length/2]  + input[i+length*3/4]
        #print '%d=%d+%d'%(i+length/2,i+length/2,i+length*3/4)
        x[i+length*3/4]     = input[i+length/2]  - input[i+length*3/4] 
        #print '%d=%d-%d'%(i+length*3/4,i+length/2,i+length*3/4)
    return x

def tf2(input):
    length = len(input)
    x=input[:]
    for i in range(length/4):
        x[i+length/4] = input[i+length/4]*w(2*i, length)
        x[i+length/2] = input[i+length/2]*w(i, length)
        x[i+length*3/4] = input[i+length*3/4]*w(3*i, length) 
    return x    

def merge_radix4_dft(input):
    length = len(input)
    x1 = []
    x2 = []
    x3 = []
    x4 = []
    x1 = radix4_dft(input[:length/4])
    x2 = radix4_dft(input[length/4:length/2])
    x3 = radix4_dft(input[length/2:length*3/4])
    x4 = radix4_dft(input[length*3/4:])
    x5 = []
    x5.extend(x1)
    x5.extend(x2)
    x5.extend(x3)
    x5.extend(x4)

    return x5
    
def radix4_dft(input):
    length = len(input)
    #print 'radix4, length is ', length
    x1 = bf1(input)
    x2 = tf1(x1)
    x3 = bf2(x2)
    if len(input)>4:
        x4 = tf2(x3)
        x5 = merge_radix4_dft(x4)
        return x5
    elif len(input)==4:
        return x3
    else:
        print 'len less than 4, error'
        os._exit(2)
        

def w(n,N):
    #print 'n is %d, N is %d\n'%(n,N)
    r = cmath.exp((-1)*1j*2*math.pi*n/N)
    #print 'r is ', r
    return r

def fft(input):
    #if len(input)!=N:
    #    print '1024 length of FFT only'
    #    os._exit(1)
    
    x1 = radix4_dft(input)
    x2 = [0]*len(input)
    for i in range(len(input)):
        x2[i] = x1[dro(i, N)]
    return x2

def ifft(input):
    x1 = numpy.conjugate(input)
    x2 = fft(x1)
    x3 = numpy.conjugate(x2)/len(input)
    return x3

f =numpy.array([
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    -1,
    0,
    0,
    -1,
    0,
    0,
    1,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0 
])

N=len(f)

t = fft(f)
g = pylab.fft(f)
tt = ifft(t)

pylab.plot(range(N), numpy.real(t),'g.',range(N),numpy.imag(t),'b.')
pylab.figure()
pylab.plot(range(N), numpy.real(tt),'g.',range(N),numpy.imag(tt),'b.')
pylab.figure()
pylab.plot(range(N), numpy.real(g),'g.',range(N),numpy.imag(g),'b.')
pylab.show()
